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We develop a lattice gas model for the nonequilibrium dynamics of microemulsions. 

Our model is based on the immiscible lattice gas of Rothman and Keller, which we refor¬ 
mulate using a microscopic, particulate description so as to permit generalisation to more 
complicated interactions, and on the prescription of Chan and Liang for introducing such 
interparticle interactions into lattice gas dynamics. We present the results of simulations to 
demonstrate that our model exhibits the correct phenomenology, and we contrast it with 
both equilibrium lattice models of microemulsions, and to other lattice gas models. 
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I. INTRODUCTION 

It is well known that oil and water do not mix, and yet, with the addition of amphiphile (or surfactant) 
to such systems, one can observe the formation of a wealth of complex structures. Eor a general review see 
Gelbart et. al. (1993). The complex structures themselves arise as a result of the particular physical and 
chemical properties of the amphiphile, the most important of these being the polar nature of the molecules, 
typically, an ionic head which is attracted to water, and a hydrocarbon tail which prefers oil. Consequently 
there is a strong energy preference for the surfactant molecules to be absorbed at, and thereby to cause the 
formation of, oil-water interfaces. The aformentioned structures occur in both binary (surfactant and water 
or oil) and ternary systems; they are in general strongly dependent on the relative quantity and nature of 
the surfactant molecules present (whether the amphiphile has an ionic or non-ionic nature, the size of the 
hydrophobic tail, and so on), and on the temperature of the system. The equilibrium properties of these 
complex structures are conveniently described by the numerous phase diagrams that have been obtained 
from experimental investigation (Kahlweit et. al 1987). The microemulsion phase itself occurs only within 
a certain region of the ternary phase space. It is of particular industrial interest because of the very low 
surface tensions that exist between the so called “middle phase” microemulsion and bulk oil and water 
phases (c/. normal bulk oil-water surface tension), see Cazabat et. al. (1982). The complex nature and wide 
industrial application of these systems make them ideal subjects for experimental, theoretical and numerical 
investigation. Most previous work in this area has been done on the equilibrium phase behaviour (Gompper 
& Schick 1995) of such self-assembling systems; comparatively little research has been devoted to their 
nonequilibrium, dynamic properties (Kawakatsu et. al 1993). Here we develop a computational technique, 
with a theoretical basis, that has the ability to simulate self-assembling amphiphilic systems under both 
equilibrium and nonequilibrium conditions. 

Lattice gases have been used as a numerical technique for modelling hydrodynamics since their introduction 
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in 1986 by Frisch, Hasslacher and Pomeau (Frisch et al. 1986) and by Wolfram (Wolfram 1986), who 
showed that one could simulate the solution of the incompressible Navier-Stokes equations using a class of 
deterministic lattice gases with discrete Boolean elements. The method was later generalised by Rothman 
and Keller (Rothman & Keller 1988) to enable the simulation of two immiscible fluids, and we use their 
model, which has since been investigated with some degree of rigour (see Sec. II, as the basic starting point 
for ours. 

The simulation of the behaviour of three-phase fluids is a complex and challenging area of computational 
science. Previously, certain lattice-gas models for three immiscible fluids (Gunstensen & Rothman 1991a) 
have been found to exhibit numerical difficulties, including a very high diffusivity of one fluid into another, 
very large interfacial fluctations and limited control over the relative strength of the surface tension coeffi¬ 
cients; such problems have encouraged a trend towards the study of simpler lattice-Boltzmann techniques 
(Gunstensen 1992). Furthermore, the complexity of microemulsion behaviour goes well beyond that of im¬ 
miscible fluids and so has been out of reach of these types of models until now. In this paper we develop 
a version of such a model that simulates the nonequilibrium, dynamical properties of amphiphilic systems, 
and in particular of microemulsions. Note that the use of the basic lattice gas model, in contrast to, for 
example, molecular dynamics simulations, gives us the ability to investigate the important late-time dynam¬ 
ics of these systems (Kawakatsu et al 1993) in a computationally tractable manner. Also, since lattice-gas 
models enable the implementation of complex boundary conditions, we can simulate very realistic interface 
formation and dynamics and investigate such systems under flow and within complex media such as porous 
rock. 

The purpose of the present paper is to define and establish the general validity of our model; we show 
that our model can reproduce certain known features of these self-assembling amphiphilic systems, we show 
that our model can reproduce. In Section we briefly review some of the work that has already been done 
on the modelling of self-assembling amphiphilic systems. Sections III and 0 contain the description and 
formulation of our model, while the results of simulations we have performed are described in Section 
These two-dimensional {2D) simulations demonstrate the ability of our model to capture the phenomenol¬ 
ogy of various self-assembling amphiphilic structures in a consistent way, including the propensity for the 
surfactant molecules to sit in thin layers at oil-water interfaces and the arrest of se paration of the immiscible 
oil-water phases when enough surfactant is present in the system. In Section we draw some conclusions 
from this work. 


II. EQUILIBRIUM MODELS OF MICROEMULSIONS 

Microemulsion models are characterized by the fundamental length and energy scales on which the system 
is described. On this basis we can categorise the three principal types of theories that have emerged, namely 
membrane models, Ginzburg-Landau theories, and microscopic models. For a review of such approaches see 
Gompper and Schick (Gompper & Schick 1995) and references contained therein. 

Membrane theories take as the basic statistical entities sheets of amphiphile, whether monolayer or bilayer, 
thereby providing a universal description for ternary and binary mixtures. The basic energy scale is that 
needed to bend the sheet, and the local free energy per unit area is written as an expansion in powers of the 
local curvature. The membrane approach is most useful for the study of the fluctations of a single membrane, 
or of those characterizing a bulk phase whose scaling behaviour is readily obtained. 

Ginzburg-Landau theories are described in terms of order parameters, which represent microscopic quan¬ 
tities averaged over extended regions of phase space. The free energy is often simply constructed from 
symmetry arguments once the choice of order parameters describing the system has been made. The ad¬ 
vantage of such theories is that they are simple enough to permit analytic analysis, to easily describe bulk 
and inter facial properties, and to calculate structure functions of the various phases. Again the mesoscopic 
length scales of the microemulsion must emerge from the theory, as they are much larger than the input 
length scales. 

The fundamental length scales of microscopic models are molecular ones. The mesoscopic lengths ob¬ 
served in the experimental analysis of such systems (scattering experiments and micrographs) must then 
be derived from the theory. All such theories to date have attempted to derive generic behaviour and so 
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the distinguishing features of particular amphiphiles and oils are largely ignored. The interactions included 
are invariably of short range, while the amphiphile molecules are treated either as scalar particles without 
internal structure or as vectors whose interaction energies depend on how they align with themselves and 
the other entities in the system. An additional simplification that is often made is to replace the continuum 
system by a discrete one, which then reduces the model to the form of a general lattice gas. Advantages 
of this lattice approach include its ability to model both bulk phase and interfacial behaviour. Further¬ 
more, the structure of interfaces, their tensions and associated wetting properties can be observed at the 
molecular level, while the equilibrium structure of the microemulsion itself can be examined by lattice-based 
Monte-Carlo simulation techniques. It is worth noting that hybrid models (Kawakatsu & Kawasaki 1990) 
for the simulation of immiscible binary mixtures with surfactant molecules have been derived; these make 
use of a combination of time-dependent Ginzburg-Landau theory for the underlying binary fluid and discrete 
particles with dipolar orientation to describe the surfactant molecules. 

Since our model is based on a microscopic treatment of the lattice gas particles, we look in more detail 
at some of the various microscopic models that have been investigated. There are numerous lattice-based 
microscopic models, including the three component model of Schick and Shih (1987), where all three com¬ 
ponents are found on lattice sites and a fairly simple Hamiltonian describes the interactions in the system. 
This includes a three-particle interaction term that allows the most obvious property of amphiphiles - their 
tendency to sit between oil and water - to be modelled without having to introduce extra degrees of freedom 
associated with their directional properties. This type of Hamiltonian can be generalized further with a four- 
particle interaction term that allows the various structures forming in binary systems, including amphiphilic 
bilayers, to be modelled (Gompper & Schick 1989). 

One much studied lattice model, introduced by Widom (Widom 1984; Widom 1986; Widom & Dawson 
1986), again has three species, but in this case the molecules are confined to the bonds of a regular lattice. 
The advantage of this model is that it is exactly equivalent to a spin-4 Ising model in a magnetic field, with 
ferromagnetic nearest-neighbour, antiferromagnetic next-nearest-neighbour, and three-spin interactions. The 
model allows a straightforward representation of many of the observed microemulsion-like structures, and it 
has the correct pattern of phase equilibria, interface phenomena and interfacial tension (Dawson 1986). There 
are many more generalisations of both the Widom and the three-component models. One such example is 
the Alexander model (Stockfish & Wheeler 1988), in which the molecules are not treated symmetrically; oil 
and water particles are placed on the lattice sites, at which up to two molecules can reside, and amphiphiles 
on the bonds, which are otherwise empty. As mentioned above, other microscopic models exist which allow 
the amphiphilic molecules to have orientational degrees of freedom. These are called vector models and 
generally consist of a simple Hamiltonian for a ternary mixture on a lattice, with an additional Hamiltonian 
that governs the orientation-dependent interactions. Some of these models allow the vector amphiphile to 
have a continuous orientational degree of freedom (Gompper & Schick 1989), while others permit the vector 
amphiphile to point only towards nearest-neighbour sites (Matsen & Sullivan 1990). 

In general the phase behaviour and/or phase diagrams obtained from all these models are reasonable, 
although for the most part they bear little actual resemblence to experimentally derived ones. This is because 
for the most part they are given as functions of theoretical interaction parameters and what one would really 
like to see is transitions due to changes in the relative concentration of amphiphile or temperature. It is 
intended that our model should have the ability to do this. The microscopic approach is also particularly 
appropriate for the study of self-assembly, in that the amphiphilic molecules are free to form interfaces or 
micelles as they choose. However, it is important to note that, aside from the hybrid model, all the models 
introduced above have the ability to model only static, equilibrium phase behaviour. The dynamics of such 
amphiphilic systems are unattainable by these methods. Another element lacking is any unified approach to 
binary and ternary systems, for example, there have been few attempts to describe a system as it evolves 
from a balanced ternary one to a binary one as oil, or water, is withdrawn. The model we present in this 
paper aims to correct many of the above mentioned deficiences. 
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III. IMMISCIBLE LATTICE GASES 


Lattice gases for the simulation of the flow of immiscible Navier-Stokes fluids were first introduced by 
Rothman and Keller in 1988 (Rothman & Keller 1988). They generalised the work of FHP to include new 
degrees of freedom and collision rules that gave rise to immiscible two-phase flow. The Rothman-Keller model 
can be described by supposing that the lattice gas particles are of two or more different colours, to denote the 
two immiscible species. They then defined the local colour flux of the particles leaving a lattice site, and the 
colour field due to particles at neighbouring lattice sites, respectively. The local colour flux is the difference 
between the two outgoing colour momenta at a site, while the colour field is defined to be the direction- 
weighted sum of the signed colour density at nearest neighbours. The Rothman-Keller collision rules were 
designed such that the “work” performed by the flux against the field is minimized, subject to the constraints 
of mass and momentum conservation. This effectively results in sending coloured particles towards regions 
of like colour, hence inducing cohesion and immiscible behaviour. Since then, lattice Boltzmann versions 
of this model have been introduced (Gunstensen et. al. 1991), and much research has been carried out on 
the theory (Rothman & Zaleski 1994), computer implementation (Gunstensen 1992), and phenomenology 
(Rothman & Zaleski 1994) of these models. In this section, we show that the Rothman-Keller model can be 
derived from an individual particle description, that will allow us to more easily motivate its generalisation 
to include surfactant particles. 

We work on a L)-dimensional lattice, £, with n lattice vectors per site. We denote the lattice vectors by q, 
where i G {1,..., n}; we note that rest particles can be accommodated in this framework by a corresponding 
zero lattice vector. The state of the Rothman-Keller model for M immiscible species at time t is then 
completely specified by the quantities ni{yi,t) G {0, ...,M}, where i G {l,...,n} and x G £. We have 
ni(x, t) = a \i there is a particle of colour a G {1,..., M} with velocity Ci at position x at time t, and 
ni{'K,t) = 0 otherwise. Thus, each site can be in any one of (M + l)’^ different states. We sometimes find it 
convenient to use the alternative representation, 

(^5^) — ^Ck:,ni (x,t) 5 

where Sa,i 3 is the Kronecker delta, though we note that the nf{'K,t) G {0,1} are not all independent since 
there can be at most one particle of any colour at a particular position, velocity, and time. 

The evolution of the lattice gas for one generation takes place in two substeps. In the propagation substep, 
the particles simply move along their corresponding lattice vectors, 

ni(x + Ci, t -f At) ^ ni{:>c,t). 

This is followed by the collision substep, in which the newly arrived particles change their state in a manner 
that conserves the mass of each species. 


n 


i 


( 1 ) 


as well as the total D-dimensional momentum. 


M n 

p(x,t) = (2) 

a i 

where we have assumed for simplicity that the particles carry unit mass. 

To further specify the collision process, we partition the (M + l)’^ different states of a site into equivalence 
classes of states that have the same values for the M + D conserved quantities (M conserved masses and 
D conserved components of momentum). Appendix A contains further information on the derivation of 
these equivalence classes. Collisions are thus required to take a state 5 to another state s' within the same 
equivalence class. If the equivalence class is of size one, as it would be if there were only one incoming 
particle, this uniquely specifies the collision process. For the more usual case in which there are many 
possible outgoing states, we must specify how to choose a single outcome. 
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We first consider the case of only two immiscible species (M = 2); we denote their colours by a = B or 
“blue” for water, and a = R or “red” for oil. In order to give these two phases cohesion, the Rothman- 
Keller model (Rothman & Keller 1988) favours collision outcomes that send particles of a given colour to 
neighbouring sites that are already dominated by that colour. To quantify this, we first define the colour 
charge of the particle moving in direction j at position x at time t, 

- nf{x,t), (3) 


and the total colour charge at a site. 


t) = (x, t) - nf{x, t)] 

j 

= 

We imagine that a colour charge q induces a colour potential, (j){r) = qf{r), at a distance r away from it, 
where /(r) is some function defining the type and strength of the potential. Its energy of interaction with 
another colour charge q' is then Hcc = q'(l){r) = qq'f{r). 

Since the collision part of the evolution process of the lattice gas conserves both the mass of each species 
and the total T)-dimensional momentum, the only contribution to the interaction energy will come from 
the propagation phase, where the outgoing colour charges do work in moving to their new sites. Hence, we 
consider the interaction energy between the outgoing particle with colour charge g'(x, t) at x G £, and the 
total colour charge g(x + y, t) at site x + y G £. If we make an infinitesimal virtual displacement of the first 
charge from x to x + ^x (not necessarily on a lattice site), the corresponding change in interaction energy is 

(5ilcc(x,t) = q'i{x,t)q{x + y,t)[f{\y - (5x|) - /(|y|)]. 

Taylor-expanding in ^x, this becomes 

SHcc{x, t) = g'(x, t)q{x + y, t)fi{y)y ■ Sx, 
where y = |y|, and where we have defined 



(^ being a positive integer or zero), giving a concise form for the derivatives of the function f(y). Since the 
outgoing particle with colour charge q[ will move in the direction c^, we let ^x —> q. We then sum over all 
outgoing colour charges i at site x, and over all sites y G £ with which they might interact to get the total 
colour work, 


AHcc{xi,t) = y2Yl t)q{x + y,t)fi{y)y ■ Ci 
i yec 



= J(x, t) • E(x, t)At, 


Ji{y)y<i{^ + y,t) At 


where we have defined the colour flux of an outgoing state. 




( 5 ) 


and the colour field. 
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( 6 ) 


E(x,t) = ^/i(y)yg(x + y,t). 
yec 


The sum over y extends over the range of the colour interaction, and will be discussed at length in Sec. IV E| . 
The above expression for the colour work is identical to that used by Rothman and Keller if the stencil 
function fi{y) is chose n to select only nearest neighbour sites. The stencil function will be discussed in more 
detail in Section IVE. For now, we note that we were able to derive the Rothman-Keller model from the 


simple assumption that the particles interact with a potential function 

Note that we have intentionally used notation that highlights the analogy of the Rothman-Keller model 
to electrostatics. Specifically, the colour charge, potential, field and flux are analogous to the electrostatic 
charge, potential, field, and the current density, respectively. The colour work can then be imagined as a 
sort of “resistive heating” resulting from moving the current of outgoing colour charges against the colour 
field in one timestep. 

The Rothman-Keller prescription is then to choose the outgoing state that minimizes Ai^cc, since this 
is the one that most effectively sends particles up the gradient of colour and thereby induces cohesion. If 
multiple outgoing states yield the same minimal value of Ai^cc? then the outcome is chosen randomly from 
among them. Rothman and Keller observed that this prescription yields phase separation when the total 
particle density on the lattice exceeds a critical value. Defining the reduced density of a lattice gas as the 
average proportion of underlying vectors at an individual lattice site that will contain a particle of some kind, 
Rothman and Keller found that the critical minimum value of reduced density needed for phase separation 
was approximately 0.2. 

This prescription was subsequently generalized by Chan & Liang (1990) who noted that the colour work, 
AHcc, can be thought of as a Hamiltonian function, H{s, s'), of the incoming and outgoing states; and they 
argued that this Hamiltonian should then be used to construct Boltzmann weights for choosing the outgoing 
state. Specifically, if C{s) denotes the equivalence class of states with the same conserved quantities as state 
s, then one can define a partition function for each equivalence class. 






( 7 ) 


s'ec(s) 


where 1//5 is a temperature-like parameter. We then define the collisional state-transition probabilities, 

( _ - _p f3H(^s,s ) if ^ { s) 

A(s^s') = i ^ ^ ’ 

f 0 otherwise 

These transition probabilities are normalized, so that 


but they generally do not obey the condition of semi-detailed balance. 


^ s') = 1, 

s 

except as ^ 0. Note that this reduces to the Rothman-Keller model as ^ oo, and to the FHP model 
with no interactions as ^ 0. Indeed, Chan and Liang observed that the phase transition to immiscibility 
occurs for a critical value of l3 as well as of density. More generally, the Chan-Liang prescription is very 
useful for constructing lattice particle simulations that include an interaction Hamiltonian iL in a manner 
consistent with the conservation of mass and momentum. 


IV. A LATTICE GAS MODEL OF MICROEMULSIONS 

To model microemulsions, we would like to introduce surfactant molecules in the framework of the 
Rothman-Keller lattice gas. To do this, we first introduce a third species index, say S, to represent the 
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presence or absence of a surfactant molecule. Thus, since M = 3, each site can be in any one of 4^^ states, 
and can therefore be conveniently represented by 2n bits. 

Surfactant molecules do not contribute to the colour flux and colour field in the same manner as ordinary 
coloured particles. Real surfactant molecules generally consist of a hydrophilic (often ionic) portion attached 
to a hydrophobic (hydrocarbon) portion. Thus, to pursue the electrostatic analogy mentioned in Sec. IH, 
surfactant particles are best imagined as colour dipoles. 

As in electrostatics, we shall model a colour dipole as a pair of equal and opposite colour charges, 
separated by a fixed displacement, a, in the limit as a ^ 0, g' ^ oo, and qa. a where a is the colour 
dipole vector. Thus, we characterize the surfactant molecule at position x moving in direction i by a colour 
dipole vector cr^(x, t). Note that the value of cr^(x, t) is zero unless nf (x, t) = 1. The total dipolar vector at 
a site is then denoted by 


n 

= y]crj(x,i). (8) 

i 

It will be necessary to take scalar products of colour dipole vectors with other vector and tensor quantities. 
Beyond this, however, we leave their precise representation unspecified for the moment. 


A. The Colour / Dipolar Field Interaction 

We first consider the work done by a colour charge, Q, moving in the field of a fixed colour dipole, cr, at 
a displacement y, as shown in Fig. 


X 


Q 


y 


-q 



+q 


FIG. 1. Model of Colour-Dipole Interaction 
The energy of the static interaction is 


Had = qQ [+/(|y + a|) - /(|y|)]. 

In the limit as a ^ 0 and g ^ oo, so that qa-^ a^ this becomes 

Hcd = -fi{y)Qcr • y. 

Since the mass and momentum of particles with colour charge g(x, t) at x G £ are the same before and after 
the collision part of the time evolution process, that is colour charge is conserved by the dynamics, the only 
local energy change that we need to incorporate into the Hamiltonian due to the above interaction comes 
about as a result of the outgoing particles doing work in moving to their new sites. To compute this energy 
change, we note that the interaction energy between an outgoing particle with colour charge q'- (x, t) at x G £ 
and the total dipolar vector cr(x + y, t) at site x + y G £ is given by 

Hcd = -fi{y)q'i{x,t)(T{x + y,t) -y 

and that if we make an infinitesimal virtual displacement of the particle from x to x + (5x (not necessarily 
on a lattice site), the corresponding change in the interaction energy, i^cd(y ~ ^x) — Hcd{y) is 
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(5iJcd(x, t) = -g'(x, i)o-(x + y,t) ■ [f2{y)yy - ■ Sx, 

where f£{y) is defined by Eq. i)- Since the outgoing particle with colour charge q'- (x, t) will move in direction 
Ci, we let Sx ^ Ci. We then sum over all outgoing colour charges i at site x, and over all sites y G £ with 
which they might interact to get the total work, Ai^cd- We find 

Ai7cd(x, t) = J(x, t) • P(x, t)At, 

where we have defined the dipolar field vector, 

P(x,i) = -'^ [f2{y)yy - fi{y)i] ■(T{x + y,t), (9) 

and where 1 denotes the rank-two unit tensor. 

Thus, we see that the effect of the dipoles at neighbouring sites is to augment the colour field felt by a 
colour charge by the amount P. The total work done by the outgoing particles with colour charge is then 

AH,, + AH,^ = J • (E + P) At. 

We must now compute the work done by the outgoing dipoles. 


B. The Dipole / Colour Field Interaction 

Next, we consider the work done by a colour dipole vector, tr, moving in the field of a fixed colour charge, 
Q, at a displacement y as shown in Fig. ||. 



Q 


FIG. 2. Model of Dipole-Colour Interaction 
The energy of the static interaction is 

= qQ [+/(|y - a|) - /(|y|)]. 

In the limit as a ^ 0 and g ^ oo, so that ga ^ tr, this becomes 

Hdc = +h{y)Q(T ■ y. 

In contrast to the case of the outgoing colour charge, the moments of the colour dipole vectors are not 
necessarily conserved by the collision phase of the update, and so the total local energy change which must 
be represented here has two parts to it. The first part, at zeroth-oidei, is due to the possibility of the 
postcollisional configuration having a different energy from the precollisional configuration, and the second, 
at first-oidei^ is an energy change resulting from the outgoing colour dipoles doing work in moving to their 
new sites. Now, note that the interaction energy between an outgoing dipole with dipolar vector cr'(x, t) at 
X G £, and the total colour charge, g(x + y, t), at site x y G £ is given by 


iJdc(x,i) = +fi{y)q{x + y,t)(Tl{x,t) ■ y. 


(10) 



and that, if we make an infinitesimal virtual displacement of the dipole from x to x + (not necessarily on 
a lattice site), the corresponding change in the interaction energy, i^dc(y — — ^dc(y )5 is 

(5iJdc(x, t) = g(x + y, t) [/ 2 ( 2 /)(t'(x, t) -yy - t)] ■ Sx, (11) 

where fi{y) is defined by Eq. ( 0 ). The zeroth order part of the total change in the interaction energy is then 
obtained directly from Eq. ( |1Q[) and, letting (5x —> c^, we can get the first order part from Eq. ([ll|). Hence 
the total change in the interaction energy, is 

SxHdc = fiiy)q{x + y,t)ai{x,t) ■y + q{x + y,t) [/ 2 ( 2 /)cr'(x,t) • yy -/i(y)<7'(x,i)] - a. 

We then sum over all outgoing colour dipoles i at site x, and over all sites y G £ with which they might 
interact to get the total dipolar color work^ Ai^dc- We find, using Eq. (^), 

Ai^dc(x, t) = t7'(x, t) ■ E(x, t) + J(x, t) : f (x, t)At, 

where the “double-dot” notation is short for Tt{J' • f), and we have defined the total outgoing dipole vector 

n 

(T\x,t) = ^(T'(x,i), 

i 

the dipolar flux tensor 

n 

J{x,t) = J2^ai{x,t), (12) 

and the colour field gradient tensor 

+ [f 2 {y)yy - fiivU], ( 13 ) 

ye£ 

and where 1 denotes the rank-two unit tensor. 

Thus, just as the change in colour interaction energy AHcc due to an outgoing configuration of colour 
charges was modelled by Rothman and Keller as the dot product of a vector flux and a vector field, we find 
that the analogous change Ai^dc due to an outgoing configuration of colour dipoles can be modelled by the 
double dot product of a tensor flux and a tensor field, together with the addition of a zeroth order term 
which depends on the colour field vector itself. To complete the model we need to look at the interaction 
energy between two dipoles. 


C. The Dipole / Dipole Interaction 

We finally consider the work done by a colour dipole, tri, moving in the field of a fixed colour dipole, t72, 
at a displacement y. To compute the static interaction energy between these two dipoles, we return to our 
fiducial model of a dipole as a pair of opposite charges, ±g, separated by a fixed displacement, a. Two such 
dipoles, separated by a distance vector y, and distinguished by subscripts 1 and 2, are shown in Eig. ||. 


+Qi 



-Qi 


y 


-q2 



FIG. 3. Model of Dipole-Dipole Interaction 
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The static interaction energy between the dipoles is then 

Hdd = qiq2 [-/(|y- ail) - /(|y+ a2|) + /(|y + a2 - ai|) + /(|y|)]. 

In the limit as a.j oo and qj ^0, so that ^ cr^, this becomes 

i^dd = -f2{y){(Ti ■ y)(cr2 ■ y) + fi{y)(Ti ■ 0-2. 

The analysis of the total interaction energy here is as in the dipole-colour field case above, since again we 
have the presence of both zeroth and first order terms to consider. Similarly, note that the interaction energy 
between an outgoing dipole with dipolar vector cr-(x, t) at x G £, and the total dipolar vector, cr(x + y, t), 
at site X + y G £ is given by 

iJdd(x, t) = -/2(2/)(cr'(x, t) ■ y)(cr(x + y,i) • y) + /i(2/)cr'(x, t) ■ (t(x + y, t). ( 14 ) 

and that if we make an infinitesimal virtual displacement of the first dipole from x to x + ( 5 x (not necessarily 

on a lattice site), the corresponding change in the interaction energy, i^dd(y “ — i7dd(y): is 

( 5 iJdd(x, t) = f2{y) {[o-'(x, t) ■ cr(x + y, t)] y + [(t(x + y,t)- y] cr'(x, t) 

+ [o-i{^,t)-y]cr{x + y,t)}-Sx- fsiy) [cr(x + y, 1 ) • y] [cr'^{x,t) • y] y • ( 5 x, ( 15 ) 

where the fe{y) are given by Eq. (^). Consequently we can use Eq. ( p^ to obtain the zeroth order energy 

change, while, since the outgoing dipole with dipolar vector ^iii iiiove in direction c^, letting ^x ^ Ci 

enables us to obtain the first order part directly from Eq. ([iq). The total change in the interaction energy, 
is then given by 

SrHdd = -f2iy){(Ti{x,t) ■y)(cr(x + y,t) ■ y) +/i(2/)cr'(x, t) ■(T{x + y,t) 

+/2(y) {[cr'(x,t) ■(T(x + y,t)]y+ [cr(x + y,t) • y] cr'(x,t) 

+ [<^i(x, t) ■ y] (t(x + y, t)} ■ dx - fsiy) [cr(x + y,t)- y] [cr'(x, t) • y] y ■ 6 x. 

We then sum over all outgoing colour dipoles i at site x, and over all sites y G £ with which they might 
interact to get the total interdipolar colour work^ Ai^dd- After some manipulation, and making use of Eq. (|y), 
we find that the result is 


Ai^dd(x, t) = a'{-K, t) • P(x, t) + J(x, t) : P(x, t)At, 
where J7(x, t) is the dipolar flux tensor, and we have defined the dipolar field gradient tensor 

P(x, i) = - E 0 • [f3{y)yyy - f2{y)y ■ , (I 6 ) 

yec 

wherein we have in turn defined the completely symmetric and isotropic fourth-rank tensor, 

^ijkl = 


D. Expression for the Total Colour Work 

The total interaction work can now be written 

A-f^int = Ai^cc + Ai^cd + Ai^dc + Ai^dd 


J + — ) •(E + P) + J:(£ + P) 


At. 
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The work done to change the kinetic energy is then 


M n 


AFke = 


[nf(x,i) -n“(x,i)]. 


( 17 ) 


The total work done is then 


^^total — ^^int + ^^ke- ( 18 ) 

We note that the change in kinetic energy was not included by Rothman and Keller in their model; neither 
has it been used in subsequent studies of their model. On the other hand, we can think of no good reason 
to omit it. Its presence or absence will not affect the equilibrium properties of the model, but may impact 
its nonequilibrium properties. 


E. Stencils 


We would like to use the total work, Ai^totai as the Hamiltonian in an implementation of the Chan-Liang 
procedure. To compute this, we will first need the colour flux and the dipolar flux tensor. These are easily 
computed directly from their definitions, Eqs. (i and respectively. Their computation involves only 
quantities that are local to the site x. 

We will also need the colour fields the colour field gradient tensor, and the dipolar field gradient tensor. 
These are given by Eqs. (H), ([l^, and respectively. Note that these expressions involve various deriva¬ 
tives of the potential function /(^), as defined in Eq. (^. In order to develop a tractable lattice model, 
however, we shall want to truncate the range of this function, and so it is unclear how we should treat 
its derivatives. In this section, we shall show how to define a single function, ^(y), which parametrizes 
the interparticle potential, and in terms of which we can derive the fields by taking linear combinations of 
quantities from some pattern of neighbouring sites. These combinations, called stencils^ involve g{y) but do 
not involve its derivatives at all. We are then free to truncate g{y) in order to develop a model that involves 
only nearest neighbour communication on the grid. 

A stencil can be specified as follows: suppose that we have a scalar field A(x). Let each site x retrieve the 
value A(x + y) from its neighbour at displacement y, multiply it by the coefficient ^(y), and sum over all 
neighbours y. The result. 


g{y)A{x + y), 

yec 

is called the scalar stencil of the scalar field A. We can also define the vector stencil^ 

Y 


and the tensor stencil^ 


Y 3(y)yy^(^ + y)^ 

ye£ 

and so on. Likewise, given a vector field B(x) we can define its scalar stencil, 

y]5(y)B(x + y), 

yec 


vector stencil. 


11 




^5(y)yB(x + y), 
yec 

and so on. We require the following isotropy properties of the coefficients g{y): 

Y, s'(y) = Go 

yec 

Y siy) y = 0 

yec 

Xifi'(y) yy = G2I 

yec 

Y yyy = ^ 

ye£ 

Y yyyy = ^4^^. (19) 

ye£ 


As a concrete example, consider a regular Bravais lattice (|ci| = c), such that all tensors up to the fourth 
order formed by the lattice vectors are isotropic, 


By choosing the stencil coefficients 


^l = n 

3 

n 

=0 

3 

n 

3 

n 

CiCiQ = 0 

3 


n 

3 


-i7. 

D{D + 2) 


9{y) = 


1 if |y| = c 
0 otherwise. 


( 20 ) 


( 21 ) 


we see at once that we have Go = n, G 2 = njD, and G 4 = n/[D{D + 2)]. 

We note that lattice gases for fluid applications must have a lattice for which Eqs. (|^ are satisfied in 
order to ensure isotropy of the resulting hydrodynamics, so it is always possible to use Eq. (H) to define a 
stencil for such a lattice gas. It is possible, however, by judicious choice of the coefficients g{y) to satisfy 
isotropy conditions like Eqs. (H) to higher rank than those of the lattice vectors, Eqs. (^q|). 


F. Stencils for the Colour Field and its Gradient 

We first note that the colour field can be computed directly from its definition, Eq. (i, in a straightforward 
manner. It is a vector stencil of the colour charge, with stencil coefficients fi{y). We could likewise compute 
the colour field gradient tensor directly from its definition, Eq. (|^, but there is a better way. Consider the 
colour field evaluated at a point x + (5x (not necessarily at a lattice site). We have 
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E(x + (5x,t) = ^ fi{\y - 5x.\)(y - 5x)q{x + y,t). 
yec 


Expanding in (5x, this becomes 


E{x + 5x,t) =E(x,i) + < ^g(x + y,t) [f 2 {y)yy - fi{y)l] > -Sx 


,yec 


+1 ^ E iMy)yyy ~ My)y • ^] |: <^x(5x. 


If we now substitute (5x —^ z in this equation, and use it to take the vector stencil of the field, we get 
E fi'(z)zE(x + z, t) = G '2 ^ q{x + y) [f 2 {y)yy - /i(2/)l], 


zEjC 


yec 


where we used Eqs. (H). Comparing this with Eq. d), we see that the colour field gradient tensor is given 
by the vector stencil of the colour field. 


^ y]fi((z)zE(x + z,i). 


( 22 ) 


zEjC 


In addition to justifying the name of this tensor, this equation is a much more efficient way to compute it in 
practice. 


G. Stencils for the Dipolar Field and Dipolar Gradient 

Einally, we need to compute the dipolar field and the dipolar field gradient tensor. Once again we could 
compute these directly from their definitions, Eqs. (i) and respectively, but there is an easier way. We 
first define the scalar field. 


S{x,t)='^fi{y)y-cr{x + y,t), (23) 

yec 

using the same stencil coefficients fi{y) that we used to get the colour field. We then consider the value of 
this field at a point x + (5x (not necessarily at a lattice site). We have 

^(x + i^x,^) = /i(|y - (5x|)(y - Sx) ■ (T{x + y,t). 

yec 

Expanding in (5x, this becomes 


5'(x + (5x, t) = 5'(x, t) + 


Y + y, i) • [h{y)yy - /i( 2 /)i] 

yec 


• ^x 


If we now substitute (5x ^ 
scalar field 5'(x, t), we get 


-1^ E ■ [/3(y)yyy - h{y)y ■ | 

z in this equation, and use it to take the scalar. 


: (5x^x. 

vector, and tensor stencils of the 
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fl(z)5'(x + z,t) = GoS{x, t) + ^Y + y^t)-y [hiy)y^ -{D + 2)f2{y)] , 
zeC yec 


^fl(z)zS'(x + z,i) = G 2 ^cr(x + y,t) ■ [f 2 {y)yy - fliy)!], 

zeC yeC 


and 


Y g{z)zzS{^ + z,t) = G2S'(x, t)i + ^1 Y + y^t)-y [f3iy)y‘^ - (d + 2)/2(y)] 


zEjC 


yec 

+ G 4 ^ cr(x + y, i) • [/ 3 (y)yyy - f 2 {y)y ■ ^ 2 ] 

y€C 


respectively, where we used Eqs. (19). Comparing these with Eqs. (H) and (16), we see that we can write 


P(x,t) = y] 5 (z)z 5 (x + z,i), 


(24) 


zEjC 


and 

^ ~ '5(x + z,i), (25) 

^ zeC ^ 

which are efficient methods to compute the dipolar field vector and the dipolar field gradient tensor, respec¬ 
tively. 


H. Sampling the Outgoing State 

We have not yet specified the exact representation of the colour dipole vectors. As has been noted, it is 
necessary that there be a colour dipole vector (x, t) associated with each particle i at each site x at each 
time t, and that it be possible to take inner and outer products of this vector with other tensor quantities. 
We now further demand that the magnitude of each of these dipolar vectors be a fixed constant a when there 
is a surfactant particle present, and be zero when there is not. That is, we have 


<7i(x,i) = anf{x,t)&i{x,t), 


(26) 


where the dipolar strength a is an input parameter, and where <t^(x, t) is a unit vector in the direction of the 
dipolar orientation. This latter quantity can be parametrized hy D — 1 real-valued angles for each direction 
i = 1,..., n, for a total of n{D — 1) angles. 

We use the Chan-Liang procedure with the Hamiltonian given by the total work, Eq. ([T^). In this 
Hamiltonian, we note that the fields depend only upon the incoming states, 5, while the kinetic energy and 
the fluxes depend upon the outgoing state, s'. More specifically, the kinetic energy and the colour flux 
depend only upon the outgoing occupation numbers, nf'(x, t), while the dipolar flux depends as well upon 
the outgoing dipolar orientations, d'^(x, t). 

With this in mind, we adopt the notation s = (n, d-), where n and & denote the occupation numbers and 
dipolar orientations of state 5, respectively; likewise s' = denotes the outgoing state. We can then 

write the Hamiltonian in the following form: 


H{s,s')=Ho{s,n') 


— ■ [E(n) + P(s)] + J(s') : [^(n) + P{s)] 


At, 


where 
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Ho{s, n') = Ai7ke(n, n') + J(n') • [E(n) + P{s)] At, (27) 

and where we have suppressed space and time dependences in favour of state dependences for clarity. Using 
Eqs. H) and (1^), we can then write 

n 

At = 

i 

where we have defined the vectors 

Ai(s) = cr [(E(n) + P(s)) + (f (n) + P{s)) ■ CjAt] (28) 


— ■ [E(n) + P(s)] + Jis') : [f (n) + V{s)] 


for i = 1,..., n. 

The partition function of Eq. is then obtained by summing over all possible outgoing occupation 
numbers and integrating over all outgoing dipolar orientations, 


Z{s) = Jd&[ - J exp 


Ai{s) 


= E< 


-dHois 


n ^ 

n J do-'iexp [-/?nfd-' ■ Ai(s)] 


The probability distribution of outgoing states is then 


P{W,&') = 


Z{s) 




The reduced probability distribution of outgoing occupation numbers (without regard to dipolar orientation) 
is then 


1 , " /■ 

^ n J ■ Ai(s)] • 

Our strategy shall be to first sample from P(n') to get the outgoing occupation numbers, n', and to then 
sample the dipolar orientations from 

^ /ws ^ exp [-(inf'&i ■ Ai(s)] 

/ d&'i exp [-pnf'&'i • Ai (s)] 

for i = 1,..., n. 

Eor example, in two dimensions we can parametrize the orientations by the angles so 

d-• = X cos 6> • + y sin . 


Writing in polar form. 


Ai{s) = Ai(s) [x cos (t)i{s) + y sin (s)], 


where 


and 


(29) 


(t)i{s) = arg{(x + zy) • ([^(n) + P(s)] • q + E(n) + P(s))} , 


(30) 
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and performing the integration we see that the partition function is 


Z{s) = (27r)”^e-^^°(*’"')n^o(/3nfAi(s)), 

n' i 

where Io{z) is the modified Bessel function. The reduced probability distribution is then 

n^o(/3nf Ai(s)). (31) 

We sample the outgoing occupation numbers n' from this, and then determine the outgoing dipolar angles 
by sampling from 


[-(inf'Ai{s) cos( 6 »' - Ms))] 

for each i from 1 to n separately. Details of these sampling procedures are given in Appendix B. Note that in 
the absence of surfactant particles n' = 0, so the Bessel functions are then all unity, and the model reduces 
to Chan and Liang’s generalisation of the Rothman-Keller model (with the addition of the kinetic energy). 


I. Summary of the Algorithm 


We can now write down the full algorithm for our microemulsion lattice gas model. Consider, for example, 
a triangular grid in two dimensions {D = 2) with up to seven particles per site {n = 7) corresponding to 
the six lattice directions plus a rest particle. We use Eq. (|^ ) for g{y), from which it follows that Go = 7, 
G 2 = 3, and G 4 = 3/4. (Note that the rest particle does not contribute to G 2 and G 4 .) The algorithm is 
then: 


1 . 

2 . 

3. 

4. 

5. 

6 . 

7. 

8 . 

9. 

10 . 

11 . 


For all sites x and lattice vectors i, propagate particle ni(x, t) and its respective dipolar angle ^i(x, t) 
to site X + c^. 


Each site x computes its colour charge g(x, t) according to Eq. (||). 

Each site retrieves its neighbours’ colour charges, and computes its colour field E(x, t) using Eq. (||). 

Each site retrieves its neighbours’ colour fields, and computes its colour field gradient tensor f(x, t) 
using Eq. (p^). 

Each site x computes its total dipole vector cr(x, t) according to Eq. (||). 


Each site retrieves its neighbours’ total dipole vectors, and computes the scalar field 5'(x, t) according 
to Eq. (H). 


Each site retrieves its neighbours’ scalar fields S', and computes its dipolar field vector D(x, t) and 
dipolar field gradient P(x, t) according to Eqs. ( p4|) and (p^), respectively. 


Each site computes its vectors according to Eq. 
to Eqs. (p^) and (30). 


I), and converts them to polar form according 


Each site uses its state s to determine its equivalance class C{s)^ and enumerate the allowed outgoing 
states s'. 


For each outgoing state 5 ', we compute the colour flux J(s') according to Eq. (||). 


For each outgoing state 5 ', we compute the kinetic energy change AiLke according to Eq. (17). 
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12. For each outgoingstate s', we compute the dipole-independent part of the Hamiltonian, i7o(<s,n'), 
according to Eq. (^) . 

13. For each outgoing configuration n', we calculate the probability P(n') according to Eq. (|M]), and 
normalize these over all outgoing configurations. 

14. The final outgoing state n' is then sampled from P(n'). 

15. Given n', the final outgoing dipolar angles ^'(x,t) are sequentially sampled from Qi{0[) given by 

Eq. (H. 

16. Go to step 1. 

For simulations using our model on moderately sized lattices in 2D we have employed the algorithm as 
it is written above. For larger system sizes, however, it is possible that this algorithm may be prohibitive 
in terms of computer memory. In this case simpler versions may be worked out in which the dipolar vector 
directions are discrete rather than continuous, and/or in which a Metropolis Monte Garlo procedure is used 
to select the outgoing state. Although the simulations reported in this paper have been undertaken using 
a basic C code, the algorithm is parallelizable, using either the message-passing or data-parallel paradigms. 
We expect to use massively parallel supercomputers to achieve larger system sizes and later times in future 
simulations. 


V. SIMULATIONS 

As mentioned earlier there are certain basic properties of self-assembling amphiphilic systems that it is 
essential our model be able to reproduce. In this section, we describe some of these features and our efforts 
to reproduce them with the lattice-gas model defined above. 


A. Common Properties of Microemulsions 

In general, the addition of a small amount of surfactant to a system of oil and water will not alter the 
two-phase coexistence; the added amphiphile will partition itself between the two phases. However, if there 
is enough amphiphile present in the system to overcome the tendency of oil and water to phase separate, 
then a fluid phase can form in which the oil and water are solubilized and the surfactant tends to be ordered 
in some way. This can result in a finite concentration of oil-in-water (o/w) or water-in-oil (w/o) droplets, 
usually called micelles^ forming within the fluid, or alternatively in the formation of sheets of amphiphile 
separating coherent regions of oil and water (sometimes called lamellar phase). If the concentrations of the 
oil and water in the system are not very different then both coherent regions will span the system, and 
the fluid is said to be bieontinuous. Fluids existing in these droplet and bicontinuous phases are termed 
microemulsions. Note that they are characterized by extensive amounts of internal interface. 

Self-assembly also occurs in the two-component binary systems of water and amphiphile. The terms micelle 
and lamellae are also used in the binary case, when they refer to an often spherical cluster of surfactant 
molecules sitting within bulk water (or an inverse micelle if in bulk oil), and sheets of water separated by 
amphiphilic bilayers, respectively. Typically, at amphiphile concentrations below a critical value, called the 
eritieal mieelle eoneentration, the amphiphile molecules exist as isolated monomers in the bulk water (or oil). 
Above the critical micelle concentration, the amphiphile molecules cluster into micelles of a characteristic 
average size, consisting of a well-defined number of monomers; although there are clusters of varying sizes 
present, the distribution of sizes is sharply peaked about one value. Lamellar phases are found at still higher 
amphiphile concentrations. 
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B. Definition of Coupling Constants 


In order to incorporate the most general form of interaction energy within our model system, we introduce 
a set of coupling constants a, /i, e, in terms of which the total interaction work can then be written as 


+ C^-^dd- 


(33) 


We carry out simulations using the algorithm described in Section |IVI|; however, we have not yet defined the 
functions fi{y) that appear in Eq. (|4) - in particular we need to specify the value of fi{y) since it appears 
in all four of the terms in Ai^int- Refering back to Eq. (0) we see that if we set 



1 if ^ = 1 
0 otherwise 


(34) 


then we retrieve exactly the equation used by Rothman and Keller for their definition of the colour field. In 
the present paper, we shall employ this simple prescription and leave any more involved usage of stencils to 
a later date. Note that we also have to specify the temperature-like parameter I3~^ (Eq. ( 0 )) , the kinetic 
energy term AiJ^e (Eq- &), and the dipolar strength a (Eq. ®)- Again, as simplifications for the current 
preliminary analysis of the model, we have set p and a to unity and the kinetic energy change Ai^ke to zero 
throughout. Unless otherwise stated, all results described below have been obtained using 2D lattices of grid 
size 128 X 128 with periodic boundary conditions in both dimensions. The colour pictures showing the time 
evolution of systems are coded as follows; black is equivalent to water, green to surfactant and red to oil. 


C. Oil-Water System 

We begin by noting that when no surfactant particles are present our model reduces to Chan and Liang’s 
generalization of the Rothman-Keller model for two immiscible fluids (Rothman & Keller 1988). With no 
surfactant particles present in the system, the only term that contributes numerically to the collision process 
is aAHcc- We set a = = 1.0 and perform a simulation that has a random initial configuration with 

equal amounts of oil and water in the system. The reduced density (i.e. the proportion of lattice vectors 
at each site that contain a particle) is 0.55. In this parameter regime, the Chan-Liang model does exhibit 
phase separation with positive surface tension, as is evident from Eig. which illustrates the nonequilibrium 
behaviour of the immiscible lattice gas. In displaying the results graphically, we have used a majority rule 
to define the displayed colour of each site; if the numbers of oil and water particles at a site are equal then 
the displayed colour is selected randomly with probability 1. The behaviour shown in the figure is, indeed, 
effectively the same as that obtained by Rothman and Keller. If left to run for a large enough time the 
system would eventually reach the completely separated state of two distinct layers of fluid. 
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FIG. 4. Nonequilibrium behaviour of a two-immiscible-Buid lattice gas, consisting of snapshots of the evolving two 
dimensional system taken at the timesteps shown. 


To make a detailed comparison between the immiscible oil-water fluid behaviour shown here and later 
simulations in which we introduce surfactant to the system, we make use of circularly averaged structure 
functions. We first calculate (Roland & Grant 1989) the two-dimensional structure factor of the colour 
charge, s(k, t). 


s(k, t) = 



qavyk-y. 


where k = (fl'KjV) (mi + nj), m,n = 1,2,q(x) is the total colour charge at a site, is the average 
value of the colour charge, L is the length of the system and N = is the number of sites in the system. 
The circularly averaged structure factor, which is what we actually evaluate numerically, is given by 


S{k,t) 


Ek ^ 


(35) 


where k = 27rnlL,n = 0,1,2,and the sum is over a spherical shell defined by (n — ^) < < 

(n+^). Note that the resolution of S{k, t) depends on kc^ the cutoff frequency associated with the lattice, that 
is, for a real-space sampling interval of A the cut-off frequency is 1/2A; above this value of the frequency there 
is only spurious information carried as a result of aliasing. In our case, kc = (27r/I/)nc, where we have chosen 
Uc to be the maximum possible value, which is half the lattice size. Figure || contains the temporal evolution 
of S{k,t) for the case of two immiscible fluids. The data is averaged over five independent simulations. 
As time increases, the peak of S{k,t) shifts to smaller wave numbers and its peak height increases. This 
behaviour is indicative of coarsening and will be used as a comparison for the surfactant-based systems 
described below. 
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FIG. 5. Temporal evolution of S(k,t). Timesteps shown are, from bottom to top, t = 0, 80, 200, 400, 800, 1400, 
2000 . 

With sufficient averaging and a rescaling of variables it can be shown (Frisch et. al. 1986) that the 
Navier-Stokes equations for an incompressible fluid are obeyed within homogeneous regions of each of the 
fluids. (There are known problems with galilean invariance due to a spurious factor in front of the inertial 
term. This issue is discussed later in this paper.) It is also possible to show, as a result of the formation 
of interfaces, that physically realistic interfacial tensions exist. In terms of the basic two-species immiscible 
lattice-gas, surface tension has been studied from both a theoretical and a numerical viewpoint in a recent 
paper by Adler, d’Humieres and Rothman (Adler et. al. 1994). Using a bubble experiment as described 
in their paper, we can investigate the validity of our basic model by evaluating the surface tension in the 
immiscible fluid case, noting that the measurement of surface tension requires the presence of an interface 
between two macroscopically defined phases. We use Laplace’s law, which in two dimensions is 

P- — P 4- — — 

Mn ^out — 

where R is the radius of the bubble, is the average pressure within the bubble and Pout average 
pressure outside. The results from simulations are shown in Fig. They give good agreement with Laplace’s 
law, and a best-fitting line through the origin results in an estimate of a 0.378. This is close to the results 
given by Adler and co-workers, although we have done significantly less averaging than in their cited paper. 
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FIG. 6. Verification of Laplace’s law and estimation of surface tension 


D. Oil-Water-Surfactant Systems 


Next we seek microemulsion-like behavior in ternary mixtures. We note that the introduction of the 
surfactant will cause all four of the terms in Eq. (|3|) to contribute in some way to the interaction energies 
of the collision process. In order that these terms enable us to reproduce basic microemulsion behaviour, we 
consider briefly what the physical contributions of the dipoles ought to be: 

1. In real microemulsions, the surfactant molecules strongly prefer to sit with their hydrophilic heads in 
water and their hydrophobic tails in oil. Since the eAHac term in Eq. (|^ ) arises as a result of the 
interaction of outgoing dipoles with the surrounding colour field, it is evident that we will want to 
include a large contribution from this term as it will favour vector dipoles aligning across, and sitting 
as near as possible to, oil-water interfaces, across which the colour field gradients are at their greatest. 

2. The term juAHcd in Eq. results from the effect of surrounding dipoles on outgoing coloured parti¬ 
cles. Since this will encourage the bending of dipoles around a central colour charge it may be important 
when it comes to looking at o/w and w/o microemulsion droplet phases, but a careful balance will have 
to be struck between this and the cAHdc term to stop it from destroying the overwhelming tendency 
for surfactant to sit at oil-water interfaces, especially when we are looking to observe bicontinuous 
structures. 

3. The final term (AHaa in Eq* arises as a result of the interactions between dipoles, and will 
allow attraction or repulsion to take place depending on the sign of (. This term may be of limited 
use for modelling more general amphiphile-containing systems because, at present, our model does 
not differentiate between the relative strengths of the hydrophobic and hydrophilic interactions in 
amphiphile molecules. 

After this assessment of the relative values of these constants, including 
• the need for e to be relatively large when compared with the other constants (see point 1. above). 


• the realisation that an effective useful maximum value for e exists, due to constrai nts im posed by the 
calculation of modified Bessel functions in the collision update process (see Section IV H| ), 
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• the need to compare our results with the known basic lattice-gas behaviour, suggesting that initially 
we maintain the value of a used in the two-immiscible-fluid case, 

• some simulation work, 

we arrived at the following set of coupling constants, that should allow us to model the various experimentally 
observed microemulsion characteristics : 

a = 1.0, iJL = 0.05, e = 8.0, C = 0.5. (36) 

We use this as a ‘canonical’ set of coupling constants in the remainder of this paper. 


E. Binary phases: From monomers to micelles 

If one adds a small amount of amphiphile to water, then those surfactant molecules, while being highly 
dynamic within the bulk water phase, will remain distinct from one another and exist as monomers. Gradu¬ 
ally increasing the amount of amphiphile in the system just increases the density of these monomers until, at 
the critical micelle concentration, the monomers begin to form micelles. These give the system characteristic 
structure and should be discernible in our simulations. The micelles themselves are dynamic objects and 
are not necessarily very long lived, since individual molecules are free to detach themselves, meet with other 
monomers and/or micelles and rejoin to form new structures; the kinetics of simple micelle formation can be 
modelled on the basis of a Becker-Doring theory (Coveney & Wattis 1995). Individual micelles do not grow 
without limit; if more surfactant is added to the system it will form new micelles rather than increasing the 
size of those already present. This is because, energetically, only a certain number of amphiphile molecules 
are able to fit around one central point to produce a micelle structure (for real micelles, the characteristic 
number of monomers contained within a micelle depends on the details of the surfactant’s molecular struc¬ 
ture). Consequently we do not expect to see evidence of micelles coalescing and growing in an unbounded 
manner in our simulations. 



FIG. 7. Time evolution of binary water - surfactant phase, for the case with surfactant - water ratio 1;8. 





phase, for the case with surfactant - water ratio 1:2. 
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We performed two simulations, one with a low concentration of surfactant and the other with a concentra¬ 
tion exceeding the critical micelle concentration of the system. In both cases, the initial condition consists 
of placing water and surfactant particles on the lattice randomly. The visual results of the simulations are 
shown below. Due to the limitations of the visualization technique employed, we require further proof of 
the existence of structure and so perform a more quantitative analysis by calculating circularly averaged 
structure functions of the surfactant density (Kawakatsu et. al Q993). For consistency, the coupling con¬ 
stants used in both simulations are as defined in Eq. Figure Q shows the result of a system containing a 
surfactant-to-water ratio of 1:8, equivalent to a initial reduced density of 0.4 for water and 0.05 for surfactant. 
As before a majority rule is employed to display the type of particle present at each site. 
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FIG. 9. Temporal evolution of surfactant density S{k,t) for binary water and surfactant mixtures. 

As can be seen in Fig. |^, with low surfactant concentration there is very little aggregation of the monomers 
throughout the simulation. This is backed up by the structure factor calculations, shown in Fig. which 
indicate no structure formation during the timescale of the run. In stark contrast. Fig. ||, which has a 
surfactant-water ratio of 1:2, initial reduced density 0.4 for water and 0.2 surfactant, clearly indicates the 
formation of small, structured objects, along with the presence of some monomers. These structured objects 
are indeed micelles; they appear in the early stages of the simulation and, although being highly dynamic, 
do seem to maintain their size and shape. This analysis is confirmed by the circularly averaged structure 
factors shown in Fig. ^ where an average taken over ten independent runs is displayed. The graph clearly 
indicates the formation of structure when the surfactant-to-water ratio is 1:2, and an absence of structure 
when the ratio is 1:8. 
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FIG. 10. Characteristic wave number {k(t)) against time t. The unbroken line on the graph corresponds to the 
case with surfactant - water ratio 1;8 and the broken line to the ratio 1:2. 


The fact that the plots at times 200 and 1000 are virtually identical suggests that the characteristic sizes 
of the micelle-like structures saturate at early times during the simulation and that they do not change 
dramatically thereafter. Thus, as expected, the structures do not grow without limit. 

This is made still more clear by Fig. |l^, in which we plot the temporal evolution of the characteristic wave 
number. 


{k{t)) 


Eko kS{k,t) 


the inverse of which is a measure of the average domain size (Kawakatsu et. al 1993). 

We see that, in contrast to the case of low amphiphile concentration, we get initial growth of the surfactant 
domains which very rapidly levels off to some constant size. 


F. Ternary phases: lamellae 

We first investigate the stability of a lamellar structure, which is composed of alternating layers of oil-rich 
and water-rich phases separated by surfactant molecules. We look at the system with and without surfactant 
present in order for a critical comparison to be made. In a similar way to Kawakatsu & Kawasaki (1990) we 
set up the initial configuration of the system, resulting in layers of oil and water eight sites wide, all sites 
having a reduced density of 0.5. It is clear that if our model is exhibiting the correct behaviour, then we 
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would expect there to be a critical density of surfactant required at the oil-water interfaces in order for the 
layered structure to be stable. Consequently, we set up a simulation where there is a layer of surfactant at 
each of the oil-water interfaces that is just a single site wide, but with a reduced density on these amphiphilic 
sites equal to 0.8. The results of the simulations undertaken are shown in Fig. Fig. 

Figure p] i 


and Fig. |13| below. 

I is the pure oil-water case with a = 1.0, ensuring that the oil and water particles will want to act 
as immiscible fluids and so we expect to see phase separation evolving from the layered initial condition. 
Figure ^ has surfactant present as described above but with coefficients a = 1.0,/i = 0.0, e = 0.0, C = 0.0, 
while F ig, has a similar amount of surfactant in the system but in this case the coefficients are as defined 
by Eq. (|36|), and so the full set of interaction terms in our model are now included. Note that we let the last 
two simulations evolve to late times in order to check that we are not just observing metastable states with 
long equilibration times which might arise as a result of the particular set of initial conditions chosen. 




FIG. 11. Time evolution from initial lamellar configuration, for the case without surfactant. 


The probabilistic, dynamical nature of our model means that there are sufficient fluctuations present to 
enable oil and water particles from the initially separated layers to move locally, and in so doing come under 
the influence of the colour field gradients produced by other layers of the same type. Since a = (3 = 1.0, 
there is an inherent tendency for the oil and water to act as immiscible fluids and phase separate (c/.. Fig. |^), 
and that is exactly what the simulation. Fig. |n|, reflects. 
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FIG. 12. Time evolution from initial lamellar configuration, for the case with surfactant but only a non-zero. 



FIG. 13. Time evolution from initial lamellar configuration, for the case with surfactant and coefficients set as 
= 1.0,/i = 0.05, e = 8.0, C = 0.5. 
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The second result, displayed in Fig. |^, is included as a control. It shows that the mere presence of a third 
type of particle at the oil-water interfaces is not enough to artificially stabilise the structure of the layers. 
However it is clear that the initial breakup of the layers takes place at a later time than for the case of two 
species only (Fig. |^). This is because of the presence of a third species between the oil and water layers, 
meaning that on average the oil/water particles will have further to move before they can come under the 
influence of a colour field set up by a neighbouring layer of like colour. An interesting point to observe here 
is the forced accumulation of the third species, in thick layers at the interfaces, once the oil and water begin 
to act as immiscible fluids and phase separate. 

Finally we come to Fig. H The presence of the colour dipoles and the fact that their corresponding 
interaction terms can affect the collision outcomes means, as the simulation clearly shows, that it is now 
energetically favourable for the surfactant to remain at oil-water interfaces. In spite of the fact that the 
dipolar vectors that represent the surfactant molecules are initially assigned random angles between 0 and 
27r, so that again we have initial fluctuations in the system, it is clear that the integrity of the lamellar 
structure is maintained and appears to be completely stable, even at late times in the simulation. At the 
oil-water interfaces the amphiphile sits in thin layers (c/. Fig. |^), which compares well with the knowledge 
that in real microemulsions amphiphile tends to reside in monolayers between regions of oil and water. The 
greater the surface tension that exists at the interfaces within a system, the more that system will act to 
try and reduce the amount of interfacial area present. Since in our model the presence of surfactant at 
the oil/water interfaces results in stabilisation of the lamellar structure, and such a structure has a large 
interfacial area, it appears that the surface tension is indeed being reduced. This analysis and result provide 
initial confirmation of microemulsion-like behaviour in our model. 


G. Ternary phases: Oil-in-water (water-in-oil) and bicontinuous microemulsions 

Finally we use our model to simulate the different ternary microemulsion phases that are possible in 2D^ 
namely, the oil-in-water droplet and the bicontinuous phases. (Note that since, at present, our model treats 
oil and water molecules in symmetric fashion, the oil-in-water phase could equally well represent a water-in- 
oil phase.) In reality the two distinct microemulsion phases will form if there are the correct relative amounts 
of oil, water and surfactant in the system at a given temperature. The oil-in-water droplet phase typically 
consists of finely divided, spherical regions of oil, with stabilising monolayers of surfactant surrounding them, 
sitting in the bulk water background. If one increases the relative amount of oil in the system and there 
is sufficient amphiphile present, one will observe the formation of mutually percolating tubular regions of 
oil and water, with layers of surfactant sitting at the interfaces. In both these cases, the equilibrium state 
does not correspond to complete separation of the immiscible oil and water regions, but rather to complex 
structures with very different characteristic length scales that form as a result of the presence of amphiphile. 

In order to reproduce the oil-in-water phase, we set up a simulation with a random initial configuration 
consisting of a 3 : 1.9 : 0.7 water-to-surfactant-to-oil ratio, respectively, with an averaged reduced density 
of 0.56. To maintain consistency between our various simulations, we again use the coupling constants as 
defined in Eq. (|3^). Figure displays the result. We see the rapid formation of many oil-in-water droplets, 
whose size increases slightly, but not without limit. This is representative of the experimental microemulsion 
state and occurs because the layer of surfactant that surrounds the droplets acts to stabilise the interfaces 
and thereby inhibits the further flow of oil to the centre of any such droplet, as well as discouraging their 
coalescence. 
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FIG. 14. Time evolution of oil-in-water microemulsion phase. 


In order to quantify this result, we calculate the circularly averaged structure factor of the colour charge, 
and plot the result in Fig. |^. We observe the early-time growth of the peak height of S{k,t) as the peak 
itself shifts towards lower values of the wavelength, indicating that the droplets form and grow to some 
characteristic size. From at least timestep 800 onwards there appears to be a negligible amount of further 
growth or movement of the position of the peak, indicating that the droplets have reached their maximum 
size and will grow no more. This simple analysis confirms the ability of our model to attain a microemulsion 
droplet phase within a certain region of the overall phase diagram. 
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FIG. 15. Temporal evolution of S{k,t) for microemulsion droplet case. Timesteps shown are, from bottom to top, 
t = 0,40,80,200,800,2000. 


To demonstrate the existence of the bicontinuous regime within the model’s phase diagram, we use the 
same coupling constants as before but simply increase the relative amount of oil present in the system. Hence 
this second simulation, shown in Fig. 16, has a random initial mixture with a reduced density 0.55 and a 


3 : 2.25 : 3 oil-to-surfactant-to-water ratio. 
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FIG. 16. Time evolution of bicontinuous microemulsion phase. 


We observe the growth of an interconnected network of tubular-like regions of oil and water, separated 
by thin layers of amphiphile that ensure the stabilisation (reduction of oil-water surface tension) of the 
bicontinuous regime, together with the formation of droplets and some micelle-like objects. Note that on 
average the width of the oil and water regions grows in size up to about 500 timesteps, during which time 
the surfactant particles migrate around the various oil-water interfaces so as to spread themselves uniformly. 
Beyond this stage the system changes very little, indicating that the observed bicontinuous phase, although 
always slightly affected by the underlying lattice gas dynamics, is stable. To appreciate the significance of 
this result, the snapshots and timescale should be compared with the two-immiscible-fluid case (Fig. |^), 
the only difference between the two being the introduction of amphiphile and the accompanying interaction 
terms. 

To permit further analysis of this result, we calculate the circularly averaged structure factor of the colour 
charge, in an exactly analagous way to that for the immiscible fluid case. The result we obtain is shown in 
Fig. 0, which is an average over five independent simulations. 
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FIG. 17. Temporal evolution of S{k,t) for bicontinuous microemulsion case. Timesteps shown are, from bottom 
to top, t = 0, 40, 80, 200, 800, 2000. 

Comparing this with Fig. ^ which shows the structure factor for the immiscible oil and water case, we see 
that after about 200 timesteps the growth of the peak of the structure factor is dramatically slowed down 
by the presence of the surfactant, and also that some minimum k is reached below which the location of the 
peak no longer shifts, indicating that the underlying immiscible fluid coalescence process has been inhibited 
owing to the presence of the amphiphile. 


VI. DISCUSSION 

The results described above provide evidence of the ability of our model to reproduce the fundamental 
equilibrium microemulsion phases. Since the dynamics conserves momentum in addition to mass, we expect 
that it ought to be able to model dynamic, nonequilibrium behavior as well. In drawing this conclusion, 
however, an important caveat should be mentioned: as noted earlier, lattice gases can break galilean invari¬ 
ance, due to the presence of a preferred galilean frame of reference, namely the lattice itself. Mathematically, 
this problem manifests itself by a spurious factor multiplying the inertial term in the momentum-conserving 
Navier-Stokes equation. For a single-phase lattice gas, this factor can easily be scaled away; for compressible 
flow, or for multiphase flow with interfaces, however, the presence of this factor is problematic, and various 
techniques have been proposed to remove it. It has been shown that this can be done at the expense of 
complicating the collision rules by introducing judicious violations of semi-detailed balance (d’Humieres et. 
al. 1987), by adding many rest particles at each site (Gunstensen & Rothman 1991b), or by using multiple 
bits in each direction rather than only one (Boghosian et. al. 1995). In the present paper, we have evaded 
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this issue by focusing only on equilibrium phenomenology, or on creep flow situations for which the inertial 
term is negligible. For flow at finite Reynolds number, however, this issue must be addressed using one of 
the above-mentioned methods, and we plan to do this in future work. 

In spite of this limitation, we note that the model exhibits some interesting effects. One of these, the 
roughening of the interface, is evident in Fig. H- The surfactant, in its attempt to increase the surface area 
surrounding a given volume of oil, has caused the creation of a fractal, fluctuating interface. Lattice gases 
have already been used as a tool for studying interface fluctuations in immiscible fluids (Adler et. al. 1994), 
we now have the ability to extend these studies to include the effect of surfactant on the interface. 


VII. CONCLUSIONS 

We have developed a model for momentum-conserving simulations of the dynamics of microemulsions and 
other related self-assembling amphiphilic systems. Using an electrostatic analogy for both the colour particles 
and the amphiphilic colour dipoles we have been able to derive the various energy interaction terms, including 
that of the Rothman-Keller immiscible fluid lattice gas, from a microscopic particulate viewpoint. Using 
a single set of coupling constants, we have shown that our model exhibits the correct 2D phenomenology 
for both binary and ternary phase systems using a combination of visual and analytic techniques; various 
experimentally observed self-assembling structures form in a consistent manner as a result of adjusting the 
relative amounts of oil, water and amphiphile in the system. The presence of enough surfactant in the system 
clearly halts the phase separation of oil and water, and this is acheived without altering the coupling constant 
a from a value that produces immisible behaviour in the case of no surfactant. In achieving these results, 
we have also demonstrated for the first time that lattice gases may be used to investigate the dynamics of 
fluids with very complex interactions. 

Consequently we should be able to investigate a plethora of microemulsion-related problems, including, 
for example, roughening and interface fluctuations in microemulsions, and the behaviour of microemulsions 
under flow conditions. Such studies are fairly easily implemented using lattice gases, and should also permit 
us to observe the various microemulsion phases flowing through complex geometries such as porous media. 
Future extensions include a ‘^D version of our model, which would allow various applications in many different 
areas of science to be investigated. 
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APPENDIX A: DERIVATION OF EQUIVALENCE CLASSES 

As discussed in the main text, the partitioning of the (M + possible states, in which each lattice site 
can exist, into equivalence classes of states that have the same values for the {M ^ D) conserved quantities 
is required for the specification of the collision process. Here we show how to derive those equivalence classes 
for lattice collisions conserving both mass and momentum. We deal specifically with the case of three species 
(M = 3) on a two-dimensional lattice {D = 2) so that we conserve three masses and two components of 
momentum. 

Since we have seven lattice vectors per site (six plus the rest particle) and each of these directions is 
assigned two bits, there is a total of n = 14 bits per lattice site and hence 2^^ possible states per site. Note 
that the two bits correspond to no molecule (00), a water molecule (01), a surfactant molecule (10), and an 
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oil molecule (11). We also need a way to identify the classes. We tag them by the number of oil, water and 
surfactant particles, the x-momentum and the y-momentum. Since there is a maximum of seven particles of 
any type at a lattice site the masses can be stored using three bits for each, while the x and y components 
of the momentum will need four bits and three bits respectively as discussed below, giving a total of sixteen 
bits of class information. 

To see why the momentum components require four and three bits, respectively, we orientate the hexagonal 
lattice so that the lattice vectors are 

. _ f xcos (^) + y sin (^) for 0 < j < 5 
I 0 for j = 6. 

Defining rii = nf^ the x-momentum is given by, 

11 11 

Px = no^ -ni - -n 2 - ns - -n^ + -ns. 

Since this can be a non-integer we multiply it by two to get a unique integer identifier for 


= 2no + ni - n2 - 2n^ - n4 + ns. (Al) 

We observe that this can be negative and that we would like an unsigned integer to identify the x-momentum, 
so we make use of the fact that 

A/^ = no + ni + n2 + ns + n4 + ns + ng, 

is a conserved quantity, since it is just the total number of particles at a site. We then add twice N to 
Eq. (^) for px to get, 

2 {^N + Px) = 4no + 3ni + n2 + n4 + 3ns + 2n6, 

which we can now use as an identifier for the x-momentum. Note that its maximum value is 14, so we need 
four bits to store this quantity. 

Similarly the y-momentum is given by 

Py = ^(^1 + n 2 - n4 - ns). 

In a similar way to that just described we multiply this equation by 2/\/3 and add getting a unique 
y-momentum identifier. 


+ —j^p^ — no + 2ni + 2n2 + ns + no. 

This has a maximum value of 7, so it can be stored in only 3 bits. 

We can now loop through all states and obtain the 16-bit class identifier for each one. These class identifiers 
are then sorted so that states in the same equivalence class are next to each other. The states and class 
identifiers can subsequently be separated into two corresponding arrays. We use the 14-bit state as a key 
to look up the initial index and length of that state’s equivalence class. All three look-up tables, the class- 
pointers, class-lengths and states are then integer arrays of length 2^^ which can be precomputed and easily 
read by the main program code at run time. 


APPENDIX B: SAMPLING PROCEDURE 

We can sample the outgoing occupation numbers n' directly from Eq. (|^). This is accomplished by 
obtaining the total energy value for each outgoing state in the particular equivalence class, normalising these 
energies and then making use of a histogram probability distribution method to select the outgoing state. 
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Similarly, we sample the outgoing dipolar orientations from Eq. Obtaining the outgoing angles from 

this, however, is not straightforward; we need to sample this function but can not do it exactly as in the case 
of the outgoing occupation numbers. The probability distribution for each dipolar angle for i = 1, • • •, n, 
is given by 

exp [-2 COS 6»"] 

' 

where 2 ; = (3nf'Ai{s) and 0" = 0[ — (t)i{s). This function is that shown in Fig. 



e" 

FIG. 18. Probability Distribution Q{0" 


For low z we use a rejection algorithm to sample from this function. This is done by choosing a random 
point, (6>, ^), uniformly distributed in the rectangle of width 27r and height e^/(27r/o(2;)). If the chosen point 
is below the curve Q{0") then it is accepted and the first coordinate of the point becomes the sampled 0"] 
if the point lies above the curve it is rejected and the procedure is repeated. One problem with this method 
is that the efficiency of the required algorithm goes as the ratio of the area under the curve to the area of 
the rectangle, so that it becomes inefficient for large 2 :. It is, however, reliable up to at least z = b. For 
larger values of 2 : we need another technique. We observe that when ^ is large the distribution Q{0") is very 
peaked near 0" = tt and can hence be approximated as follows: 

exp [- 2 ; cos 6>"] 
exp [z cos {0" — tt)] 


exp 


where we have used a Taylor expansion about {0" — tt) and the asymptotic expansion of Iq{z). This gives us 


27tIo{z) 


Q{0") = 



We observe that this is a normalised Gaussian in 0", centered at tt and of width y/ljz^ and consequently it 
is straightforward to sample 0" from this using the 2D Box-Mueller algorithm. 
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